Rayleigh loops in the random-field Ising model on the Bethe lattice 
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We analyze the demagnetization properties of the random-field Ising model on the Bethe lattice 
focusing on the beahvior near the disorder induced phase transition. We derive an exact recursion 
relation for the magnetization and integrate it numerically. Our analysis shows that demagnetization 
is possible only in the continous high disorder phase, where at low field the loops are described by 
the Rayleigh law. In the low disorder phase, the saturation loop displays a discontinuity which is 
reflected by a non vanishing magnetization after a series of nested loops. In this case, at low 
fields the loops are not symmetric and the Rayleigh law does not hold. 
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I. INTRODUCTION 



A ferromagnetic material is characterized by a rema- 
nent magnetization even at zero field. In several in- 
stances, however, it is convenient to demagnetize the sam- 
ple, bringing it to a state of zero magnetization at zero 
field. In practice, this can done applying a slowly vary- 
ing AC field, decreasing its amplitude after each cycle. 
In this way, the system explores a complex energy land- 
scape, due to the interplay between structural disorder 
and interactions, until it is trapped into a low energy 
minimum. If the demagnetization process is performed 
adiabatically and thermal effects do not play an impor- 
tant role, the demagnetized state is reproducible for a 
given realization of the quenched disorder and can thus 
be used as a reference to define magnetic properties. 

The hysteresis loops at low fields, starting from the de- 
magnetized state, are usually described by the Rayleigh 
W law: when the field is cycled between ±H* , the mag- 
netization to follows to = (a + bH*)H±b((H*) 2 -H 2 )/2 1 
where the signs ± distinguish the upper and lower branch 
of the loop. Consequently the area of the loop scales with 
the peak field H* as W — 4/36(if*) 3 and the response 
to a small field change, starting from the demagnetized 
state is given by M* = a{H*) ± b(H*) 2 |,|. This law 
has been measured in a variety of materials, but a few 
papers have reported significant deviations from the sim- 
ple quadratic law but no explanation has been provided 
I- 

The current theoretical interpretation of this law is 
based on a 1942 paper by Neel M, who derived the law 
formulating the magnetization process as the dynamics 
of a point (i.e. the position of a domain wall) in a ran- 
dom potential. In this framework, the initial suscepti- 
bility a is associated to reversible motions inside one of 
the many minima of the random potential, while the hys- 
teretic coefficient b is due to irreversible jumps between 
different valleys. The main drowback of Neel theory relies 
in its purely phenomenological nature, being based on a 
zero dimensional model which does not include collective 



effects considered very important for the magnetization 
process |6[||. 

In the past few years, the zero temperature random- 
field Ising model (RFIM) has been used to describe the 
competition between quenched disorder and exchange in- 
teractions and their effect on the hysteresis loop ||. In 
three and higher dimensions, the model shows a phase 
transition between a continuous cycle for strong disorder 
and a discontinuous loop, with a macroscopic jump, at 
low disorder. The two phases are separated by a sec- 
ond order critical point, characterized by universal scal- 
ing laws [p|-pl| . A behavior of this kind is not restricted 
to the RFIM but has also been observed in other mod- 
els, with random bonds [ p2[ , random anisotropies |l3| ] or 
vectorial spins |l4| . In addition, a similar disorder in- 
duced phase transition in the hysteresis loop has been 
experimentally reported for a Co-CoO bilayers (l^] . 

The RFIM is probably the simplest model including 
disorder and exchange interactions that can be treated 
analytically. The model has been solved exactly in one 
dimension jl6 17 1 and on the Bethe lattice ||l8| - ^o||, while 
mean- field theory § and renormalization group |10|] have 
been used to analyze the transition. Recently the one di- 
mensional solution of the model, has been generalized to 
obtain the complete demagnetization process and to de- 
rive the Rayleigh loops ||l[] . The RFIM does not display 
a phase transition in one dimensions, while numerical 
simulations indicate that the transition has an impor- 
tant effect on the demagnatization process. In partic- 
ular, in the low disorder phase the discontinuity in the 
saturation curve prevents the magnetization to reach a 
demagnetized state |^l[ but this behavior has not been 
understood theoretically. It has been shown exactly that 
the RFIM displays a disorder induced phase transition on 
the Bethe lattice when the coordination number z > 4 
JLsl p0[ and can thus be used to clarify the issue. 

Here we generalize the analysis of Refs. |Ls|-pO| and H] 
to obtain exact recursion relations for the demagnetiza- 
tion process on the Bethe lattice and show that demag- 
netization is only possible in the high disorder phase. In 
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the low disorder phase the remanent magnetization after 
a series of nested loops of decreasing amplitude does not 
vanish but scales to zero as the transition is approached. 
Furthermore, in the low disorder phase the Rayleigh law 
is not obeyed and low field loops are not symmetric. The 
Rayleigh law is instead recovered in the high disorder 
phase and the Rayleigh parameters a and b behave qual- 
itatively as in d — 1, displaying a peak in the disorder. 

The paper is organized as follows: in section II we de- 
scribe the model. In section III we recall the results ob- 
tained in Refs. p^-pp[| and generalize them for a series of 
nested loops. Section IV discusses the effect of the phase 
transition on the magnetization and section V analyzes 
the Rayleigh law. A brief discussion of the perspectives 
is reported in section VI. Finally, the complete derivation 
of the recursion relations is reported in the Appendix. 

II. THE MODEL 

In this section we recall briefly a model used to describe 
hysteresis loops in magnetic materials 11]: the ferromag- 
netic Random Field Ising Model (RFIM). This model is 
characterized by the Hamiltonian 

H = - J} j SjSj - H) Sj - hjSj (1) 

where J > 0, the ^ is restricted on the pairs of near- 
est neighbors on a lattice of coordination number z, s$ is 
the Ising spin on the site i, H is a homogeneous exter- 
nal field, and hi represents a quenched random field on 
the spin Si modelling the presence of lattice defects. The 
fields {hi} are independently drawn from a symmetric 
distribution p(hi). In the following the numerical results 
are referred to a Gaussian distribution with variance R 2 . 

In this paper we study the case of a Bethe lattice with 
a generic coordination number z. In particular we are 
interested in the case z = 4 which is known to be the 
minimal case showing a disorder induced phase transi- 
tion towards a discontinuous hysteresis loop. 

In order to mimic the microscopic spin dynamics, we 
use the flipping rules used in Refs. J|-[ll],|l6| go) obtained 
from the Glauber dynamics at temperature T and with 
an external field of frequency u> taking the limit T —> 
first and then w — > 0. The basic rule of this T = 
dynamics is that the spins align with the local field: 

Si = sign(/ij, e //) , (2) 
where the effective local field felt by the spin i is 

hi, e ff =-— =J ^2 Sj + H + hi (3) 

and the sum runs over the z nearest neighbors of the site 
i. 



Note that though the model is defined through three 
external parameters J, H, R, the dynamics is determined 
by the two reduced quantities H/R and J/R only. For 
the sake of simplicity, from now on we rename these two 
ratios H and J respectively and consider R = 1. Given H 
and J, we can write the probability p m (H) that a spin i, 
with m (0 < to < z) of its neighbors up, is also up. This 
is given by the probability that /ii, e // > 0: 

/>OC 

p m (H) = P(h i>eff > 0) = / dh'p(h'). (4) 

J(z-2rn)J-H 



III. HYSTERESIS LOOPS 

In general a change in the applied field H produces a 
rearrangement of the spins, so that each spin i is sta- 
ble being aligned with its effective field fri, e //- It is im- 
portant to note that each spin flip modifies the effective 
field on the nearest neighbors and sometimes generates 
an avalanche of spin flips through the lattice. In the 
following we will consider the case of a slowly varying 
external field: its value is kept constant until the next 
metastable state is reached. Two important properties 
of the T = dynamics are (i) the Abelian property - 
the stable state after an avalanche does not depend on 
the order in which the spins flip -, and (ii) the return 
point memory - when the field is changed adiabatically 
the stable state only depends on the point where the field 
was last reversed. These two properties can be used to 
obtain exactly the shape of the hysteresis loops. We first 
recall the derivation of the saturation curve and those of 
the first minor loops and then procede with the general 
analysis of minor loops. 



A. Saturation loop 

When the external field H is cycled from — oo to +oo 
and back the magnetization describe the saturation loop. 
The key quantity describing the lower half of the satu- 
ration loop is the conditional probability Uq{H) [H de- 
fined as the probability that a given spin flips before a 
fixed nearest neighbor, conditioned to this neighbor be- 
ing down. The probability Uq(H) satisfies the following 
equation Jl8| : 

MH) = E ( Z ~ [Uo{H)] m [l - UoiH^-^p^H) . 

(5) 

It can be shown that the probability that a spin is up at 
external field H is 
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= E j ^o(ff)] ro [l - C/o(i?)] 2 -> m (i?) ■ (6) 

m=0 ^ ' 

The related magnetization is given by mi(H) = 1 — 
2p(H), which describes the lower half of the hysteresis 
loop. The upper half of the hysteresis loop can then be 
obtained by symmetry (i.e. m u (H) — mi(—H)). 



B. First minor loops 

If the external field H is raised from — oo to a finite 
value Hq (i.e. we are on the lower half of the major loop) 
and then it is reversed, the magnetization describes the 
upper half of a minor hysteresis loop. When the field is 
reversed from Hq to Hi < Hq we define the conditional 
probability D\(Hi) for a spin to be down before a fixed 
nearest neighbor, conditioned to this neighbor being up. 
The probability Di{Hi) satisfies the following equation 




£>i(ffi) = E f 1 )[Po(Ho)r{[l-Po(H»)] ! - 1 - m [l 

m=0 ^ ' 



The related probability that a spin is up at H± is 



Equation (Q) holds as long as H± is larger than Hq — 2 J 
pofl . In one dimension it has been shown |l7]] that at 
H = Hq — 2 J the upper half of minor loop merges with 
the upper half of the major loop with the same local 
slope. This proof can be extended to the case of a Bethe 
lattice as long as the saturation loop is continuous p0| . 
The case in which the saturation loop displays a discon- 
tinuity is discussed below for the case z = 4. 

C. General formula for nested loops 

The method used to find Uq and D\ can be general- 
ized to obtain a complete characterization of all minor 
loops. In particular we are interested in nested minor 
loops, since they are directly related to the demagnetiza- 
tion process of the disordered ferromagnet. Nested loops 
is defined as follows: after having reached Hi, we reverse 
again the field increasing its value up to H2 < Hq (lower 
half of the first minor loop). This process is then it- 
erated in a sequence of fields i?2n G [i?2n-i, flzn+l] an d 
H 2n +i £ [H 2n , H 2n +2\ with n > 1, where H 2n and H 2n +i 
refer to the final value of the field H in the lower half of 
the n minor loop, and for the upper half of the (n+ l) th 



-p m+1 (H )} + [D^H^Y-^^+^Hq) . 

(7) 

(8) 



minor loop respectively. The generalizations of Uq to the 
minor loops is called V 2n (H2n) while the generalization 
of Di(Hi) is called D2 n +i{H 2 n+i) ■ In what follows, we 
simply indicate U2n(H 2n ) with C/ 2 „ and D2n+i(H 2n +i) 
with D 2n +i- 

Since H 2n < #2n-2, the set of spins contributing to 
U2n will be a subset of those contributing to U~2n-2, so 
that we can write 

U 2n = U 2n -2 - Tfcn-l + V2n (9) 

where r\2n-\ represents the fraction of spins that were up 
at i?2n-2 before their fixed nearest neighbor and down at 
H2n-ii while ri2n is the fraction of the set contributing 
to r\2n-\ which flip up again at i? 2n . The explicit deriva- 
tion of r\2n-\ and r\2n is a little involved and it is thus 
discussed in the appendix. 

The magnetization at H = H2n is as usual obtained 
as m2n = rn(H2n) = 1 — 2p(7?2n), where the probability 
p{H2n) that a spin is up at i?2n is given by the probabil- 
ity p(H2n~i) that it was already up at i?2n-i summed to 
the probability to flip up when the field goes from H 2n _\ 
to H 2n : 



m—0 



p(H 1 )=p(H ) - [ Z )[UQ(HQ)} m [D 1 (H 1 )r- m [p m (H ) -pM)] 



r 



m=0 



(fl2n) =P(H 2 n-l) {U2n} m [D2n-l} Z - m [p m (H2n) - P m (H 2n -l) 



(10) 
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The generalization of D\ to nested minor loops is called 
D2n+i and analogously to U2n{H2n) is given by 

D 2 n +l = -D2n-1 ~ C2n + Qln+l ? (H) 

where is the fraction of spins that were down at i?2n-i 
before their fixed nearest neighbor and at -ff2n and 
C2n+i is the fraction of the set of spins contributing to 
(2n which flip down again at H 2n . The exact expres- 



In principle, an arbirary series of nested loop can be ob- 
tained solving the recursion relation for U-zn, -D271+1, and 
using the result to obtain the magnetization. A similar 
procedure was used in d = 1 to obtain a closed expression 
for the magnetization along the demagnetization curve. 
This is not possible for the Bcthe lattice where an explicit 
solution for the problem is not available and one should 
resort to a numerical integration. 

IV. DISORDER INDUCED PHASE TRANSITION 
AND DEMAGNETIZATION 

Previous numerical studies of the zero temperature dy- 
namics of the RFIM on a regular lattice in finite dimen- 
sion (|Jll[] have shown that in d > 3 the system shows a 
phase transition from a strong disorder phase to a weak 
disorder phase, separated by a second order critical point 
(in d = 2 the presence of a phase transition is still contro- 
versial). In the strong disorder phase the major hysteresis 
loop is continuos, whilst in the weak disorder phase the 
loop shows a macroscopic jump in the magnetization at a 
critical value of the field. Note that if we use the reduced 
parameters J/R and H/R, fixing R = 1, strong disorder 
corresponds to small values of J and the phase transition 
will be charcterized by a critical value J c of the exchange 
coupling J. On the Bethe lattice, a phase transition is 
observed for coordination numbers z > 4, while for z < 3 
one has only the strong disorder phase M|. The pres- 
ence of the phase transition has strong implications on 
the possibility of demagnetizing the system. In particu- 
lar in the weak disorder phase phase it is not possible to 
demagnetize the system through an oscillating external 
field with decreasing amplitude. 

Before we specialize to the case z = 4, let us note 
that for H = J one has p z -i- m (J) = 1 — p m (J). This 
allows, after a little algebra, to show that for H = J, 
Uq(J) = 1/2 is a solution for any J, any z and any ran- 
dom field distribution. For z = 4 Eq. (||) is a cubic 
equation in Uq, with coefficients depending on H and J 
through the pi. In order to find the critical point it is 
enough to find the value of H and J corresponding to 



sion for the fractions £2™ and C271+1 are reported in the 
appendix. 

The magnetization at H — H 2n +i is rn 2n = m(H2 n ) = 
1 — 2p(H 2n +i), where p(H 2n +i) ^ s the probability for a 
spin to be up at H 2n +i- This probability can be writ- 
ten as the analogous probability at H = H 2n minus the 
probability to flip down between H 2n and H 2n+ \: 



(12) 



a triple solution of the equation. Implementing this re- 
quirement for any simmetric density function of the dis- 
order, one finds H c = J c , where J c satisfies the equation 
Po{J) +Pi{J) — 1/3. For a Gaussian distribution of the 
disorder, this translates in the following implicit equation 
for J c : 

Erf(J c ) + Erf(3J c ) = 1/3 (13) 

resulting in J c = 0.56140099587319..., in accordance with 
the result quoted Ref. [ p0[ . 

Above the transition (J > J c or weak disorder phase) 
the hysteresis loop becomes discontinuous Jig ]. In fact 
at J = J c the and H = J c the susceptibility dm/dH 
diverges, and for J > J c one observes a discontinuity 
whith a spinodal singularity. At this point one can mea- 
sure a gap Am in the magnetization. It is easy to show, 
through an expansion of Eq. (j^) to the lowest order in 
J — J c around J c , that for J — » it is 

Am - (J - J c f (14) 

with (3 — 1/2 as in the mean field case. The analyt- 
ical derivation of this result can be easily sketched as 
follows: first of all the three solutions of Eq. (||) at 
H = J can be found explicitly. They areU^ = 1/2 
and^ fc) ' (c) = 1/2(1±^(1 - 3(po +pi))/(l - 3pi + Po )), 
where, for H = J, (1 — 3(po + Pi)) > only for J > J c 
and = at J = J c . The gap can be measured by 
AU =| Uq — Uq C ^ I, and it is simple to show that 
AU ^ ( J — Jc) 1 / 2 , which gives Eq. ( p^ ) considering that 
to the lowest order Am ~ AU. Actually, the position of 
the gap in the mgnetization of the major loop of the hys- 
teresis cicle is located at H g > J, which could give correc- 
tions to the previous result. However, one can show that 
(H g — J) ~ (J — J c ), implying corrections to the previous 
value of the gap of the same order in (J — J c ). Then this 
correction does not alter the found scaling behavior. 

We recall that usually it is possible to demagnetize a 
material by applying a slowly oscillating external field ap- 
propriately chosen. In practice, this corresponds to a se- 



p(H 2n+1 ) = p(H 2n ) - J2 ( Z )[U2n} m [D 2n+1 } z - m [p m (H 2n ) - p m (H 2n+1 )} 
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ries of nested loops, starting from the completely magne- 
tized situation (e.g. m = —1 and h — > — oo) and then ap- 
plying in succession the fields Hq — J, Hi = —Ho(l — e), 
H 2n+ i = —H2 n (l — e) in the limits e — > + and 
n — » +00 pi] ]. This sequence, for J < J c , leads to a 
completely demagnetized state. This is due to two fun- 
damental properties of the strong disordered phase: (i) 
the "return point memory" property which has been de- 
fined in the previous paragraph; (ii) the fact that, if we 
are on a certain point of the saturation loop (e.g. at 
H = Jo on the lower half) and invert the field to Hi 
the system meets the other half on the saturation loop 
if Hi = Ho — 2J. This implies that in order that the 
first minor loop is symmetric with respect to the origin 
H = and m = without touching the saturation curve 
we have to start from Hq = J on the lower half of the 
saturation loop (or equivalently from Hq = —J on the 
upper half). 

At J > J c the demagnetization process is no more 
possible, because the discontinuity prevents minor loops 
to be symmetric with respect to the origin of the axes. 
In fact there is now an inaccessible region of the plane 
(if, m) around the origin H]. However, the field suc- 
cession described above still provides a well defined pro- 
cedure (apart the broken simmetry m — > — m and H — > 
—H) to minimal possible residual magnetization, that we 
denote m,oo, at H = 0. Following this procedure, by nu- 
merically integrating Eqs. (p|,pd|) with e = 10~ 3 , we find 
that moo displays the same scaling behavior of the gap 
in the saturation loop as J approaches J c from above: 
TOqo oc ( J — Jc) 1 / 2 (see fig(^|)). Notice that in numerical 
simulations in the d = 3 RFIM Q moo was found to 
coincide, in the low disorde phase, with the saturation 
magnetization too(0). 

V. RAYLEIGH LAW 

The Rayleigh law describes the hysteresis behavior at 
low field in a vast class of materials. Exact values for 
the Rayleigh parameters have been obtained for the one 
dimensional RFIM where the initial susceptibilty a 
and the hysteretic coefficient b both display a peak in the 
disorder R. A similar behavior is observed in simulation 
for d = 2, 3 but only in the high disorder phase, while in 
the low disorder phase a and b are not defined. 

In the case of the Bethe lattice, we could not obtain 
an explicit expression of the Rayleigh parameters even 
in the high disorder phase. We thus resort to numeri- 
cal integration and analyze the demagnetization curves 
close to H = 0. We estimate the susceptibility a and 
the hysteretic coefficient b using a linear fit of m 2n / H 2n 
vs H-2n- According to the Rayleigh law for n — > 00 we 
have m 2n /H 2n = a + bH 2n , and similarly for negative 
fields m 2n +i/H 2n+ i = o, — bH 2n+ i. The values of a and 
b as a function of the exchange coupling J are shown in 



Figs H and 0. When plotted as a function of the disorder 
R, a and b show a peak in the high disorder phase, in 
agreement which the results on Euclidean lattices. 

In the low disorder phase the demagnetization curve is 
not symmetric with respect to H — and consequently 
the Rayleigh law does not hold. In particular, we can 
define two values for the coefficient b: 

m 2n /H 2n = a + b + H 2n and 

m 2 n+i/H 2 n+i = a - b~H 2n+ i. (15) 

The values of a, 6 + , and b~ as a function of the exchange 
coefficient J are shown in fig (|],[|); also the difference 
A6 = b + — b~ is shown in fig (||), Again, as J — » J^ , Ab 
approaches as (J — Jc) 1 / 2 (see fig@). 

VI. DISCUSSION 

In conclusion, the present analysis allows to clarify the 
role of a disorder induced phase transition on the demag- 
netization properties of a ferromagnet. In particular, we 
find that in the low disorder phase the jump in the sat- 
uration curve gives rise to an inaccessible region in the 
(m,H) plane close to if = and m = 0. Even after an 
infinitesimally fine series of nested loop the final magneti- 
zation Woo does not vanish. Approaching the transition, 
however, moo scales to zero with an exponent 1/2, which 
is the same as the one controlling the size jump in the 
saturation curve. While this value could be an artifact of 
the Bethe approximation, the impossibility to demagne- 
tize the system in the low disordered phase has already 
been observed in three dimensional numerical simulations 
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APPENDIX 

Here we derive the expressions for rj 2n ^i, r] 2n , C,2n and 
C2n+i appearing in Eqs. (|J) and (|TT|) . We first note that 
these quantities can be defined in a recursive way. In par- 
ticular, r\ 2n -i represents the fraction of the set of spins 
which, at H 2n ^ 2l contribute to U 2n - 2 , but are down at 
H 2n ~i- To obtain the weight associated to this fraction, 
consider a spin i with a given neighbor j kept down (i.e. 
the spin j is conditioning the probabilities). When the 
spin i flips down at H 2n _i, apart from j all the other 
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neighbors can be either up or down. Consider for in- 
stance the case in which m of these neigbors are up and 
z — 1 — m down. Under these conditions, the spin i flips 
down if its effective field is positive at H2 n -2 and nega- 
tive at Hin-x < H2n-2- The associated contribution to 
f]2n-\ is then given by 

[U2n^2] m [D2n-l] Z - 1 - m [p m {H2n-2) ~ Pm(H 2 n-l)} (16) 



where [C/2n-2] m [£ , 2n-i] z_1_m is the probability that, at 
the moment at which the spin i flips down, m given neigh- 
bors are up and z — 1 — m (other than j) are down, 
and \p m (H 2n - 2 ) —Pm(H 2n -i)] is the probability that the 
spin i, having m up neighbors, is up at H 2n — 2 but not 
at H2n-x- To obtain r^n—ii w ^ have first to multiply 
Eq. ( |ilf ) by a combinatorial factor ( z ^ 1 ), taking into ac- 
count all the equivalent choices of the site j, and then 
sum over m from to z — 1 . The results reads: 



V2n-1 = { ' )\U2n-2Y n [D2n-x] Z - l - m \p m {H2n-2) - Pm (H 2n ^)]. (17) 



An analogous procedure can be implemented to derive 

the fraction of the set of spins contributing to V2n ^ ^ ong discussed abov6j we obtain 
which flip back up at H2 n . Using a derivation similar to . 



m— 



z-l 



[u 2 



l [p m (i?2n) -Pm(H 2n -l)]- 



(18) 



The quantities C,2n and C2n+i can be obtained proceeding 
as in the evaluation of r]2n-\ and 77271, noticing that in this 
case the fixed neighbor j (conditioning the probabilities) 
has to be kept kept up. Moreover the spin i must flip 



from down to up at H2m in order to contribute to £2™, 
and then flip back down at H2 n +i in order to contribute 
also to C2n+i- The final results reads 



( 2 n = El=0 ( Z m 1 )[U2n] m [D2n-l] Z 1 ™ [Pm+1 (H 2n ) - Pm +1 (H 2n - 1 )] 
C2n+1 = El=0 ( Z m 1 )[^2n] m [-D2n+l] Z_1_m [Pm+l(^2n) - Pm+1 (i?2n+l )] 



(19) 
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FIG. 1. Final magnetization m^o as a function of the ex- 
change coefficient J. 
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FIG. 3. Susceptibility a as a function of the exchange co- 
efficient J. 
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FIG. 4. Hysteretic coefficient b as a function of the ex- 
change coefficient J. 
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FIG. 2. The Rayleigh parameters (solid line is a, dashed 
line is b) as a function of the disorder width R for 
R > R c = \ j J c . 
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FIG. 5. The deviations from the Rayleigh law are measured FIG - 6 - Final magnetization and Ab as a function of 

by Ab = b + — b~ . J — J c . The solid line has a slope 1/2. 
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